R1 Indexes Cortical Myeloarchitecture

Read in Data for Myeloarchitecture Characterization

Final study sample

participants <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/sample_info/7T_MP2RAGE_finalsample_demographics.csv") #215 scans from 140 subjects. csv generated by /sample_construction/finalsample_7Tmyelin.Rmd

Glasser regions list

glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")

Depth-dependent R1 in HCP-MMP regions

#R1 by region matrices at 7 cortical depths
myelin.glasser.7T <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/depthR1_glasseratlas_finalsample.RDS") #generated by /surface_metrics/surface_measures/extract_depthdependent_R1.R
ordered_depths <- c("depth_7", "depth_6", "depth_5", "depth_4", "depth_3", "depth_2", "depth_1")
depth_colorbar <- c("#004A38FF", "#006458FF", "#0093B0FF", "#8AABD5FF", "#C0BFE3FF", "#E0D5ECFF", "#F2E7F4FF")

Regional average R1 in HCP-MMP regions

#Calculate across-depths average regional R1
regionalaverage.myelin <-  do.call(rbind, myelin.glasser.7T) #R1 at all depths
regionalaverage.myelin <- regionalaverage.myelin %>% select(contains("ROI"))
regionalaverage.myelin <- colMeans(regionalaverage.myelin) %>% as.data.frame()  %>% set_names("mean.R1")
regionalaverage.myelin <- regionalaverage.myelin %>% mutate(orig_parcelname = row.names(regionalaverage.myelin))

Myelin Maps

#Myelin Water Fraction Imaging
MWF.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/MyelinWaterFraction_Liu2019/source-Liu2019_desc-MWFmap_space-fsaverage_den-164k_atlas-glasser360.pscalar.nii")
MWF <- data.frame(MWF = MWF.cifti$data, orig_parcelname = names(MWF.cifti$Parcel))

#Myelin Basic Protein Gene Expression
MBP.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/AHBA_magicc/expression_maps/source-magicc_desc-MBPexpression_space-fsLR_den-32k.pscalar.nii")
MBP <- data.frame(MBP = MBP.cifti$data, orig_parcelname = names(MBP.cifti$Parcel))
MBP <- rbind(MBP, MBP) #reflect gene expression across hemispheres
MBP$orig_parcelname[1:180] <- glasser.regions$orig_parcelname[1:180]

#Sensorimotor-Association Axis
SAaxis.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii")
SAaxis <- data.frame(SA.axis = rank(SAaxis.cifti$data), orig_parcelname = names(SAaxis.cifti$Parcel))
df.list <- list(regionalaverage.myelin, MWF, MBP, SAaxis) #dfs to merge
myelin.maps <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
myelin.maps <- merge(myelin.maps, glasser.regions, by = c("orig_parcelname"), sort = F)
myelin.maps <- myelin.maps[!(myelin.maps$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

R1 Indexes Cortical Myelin Content Across Cortical Regions

myelin.colorbar <- c("#edf6ff", "#c2d4ed", "#809dc4", "#345c91", "#0e4669")

Whole-brain Map of R1

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), , colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.54, 0.58), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/R1_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

R1 Correlation with Myelin Water Fraction Imaging

MWF Map

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MWF), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.02, 0.05), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

Correlation

cor.test(myelin.maps$mean.R1, myelin.maps$MWF, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  myelin.maps$mean.R1 and myelin.maps$MWF
## S = 2261576, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6422757

Correlation Plot

myelin.maps %>% filter(MWF < 0.08) %>% #remove 2 MWF major outliers in PHA1
ggplot(., aes(x = MWF, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nMyelin Water Fraction", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

R1 Correlation with Myelin Basic Protein Gene Expression

MBP Map

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MBP), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(-0.65, 1), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

Correlation

cor.test(myelin.maps$mean.R1, myelin.maps$MBP, method = c("spearman"))
## Warning in cor.test.default(myelin.maps$mean.R1, myelin.maps$MBP, method =
## c("spearman")): Cannot compute exact p-value with ties
## 
##  Spearman's rank correlation rho
## 
## data:  myelin.maps$mean.R1 and myelin.maps$MBP
## S = 3214088, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4916123

Correlation Plot

myelin.maps %>% 
ggplot(., aes(x = MBP, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nMyelin Basic Protein Expression (z)", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

R1 Correlation with the Sensorimotor-Association Axis

S-A axis

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = SA.axis), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, trans = 'reverse', na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

Correlation

cor.test(myelin.maps$mean.R1, myelin.maps$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  myelin.maps$mean.R1 and myelin.maps$SA.axis
## S = 10280562, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.6261257

Correlation Plot

myelin.maps %>%
ggplot(., aes(x = SA.axis, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nSensorimotor-Association Axis", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

R1 Indexes Cortical Myelin Content Across Cortical Depths

#Compute mean R1 in each region at each cortical depth
regional_means <- function(input.df){
  meanmap <- input.df %>% select(contains("ROI")) %>% colMeans() %>% as.data.frame() %>% set_names("mean.R1")
  meanmap$orig_parcelname <- row.names(meanmap)
  meanmap <- merge(meanmap, glasser.regions, by = c("orig_parcelname"), sort = F)
  return(meanmap)
}

regional.myelin.alldepths <- lapply(myelin.glasser.7T, function(depth) {
  regional_means(depth)
})

Regional R1 Distributions by Depth

regional.myelin.alldepths <-  do.call(rbind, regional.myelin.alldepths) #regional mean R1 at all depths
regional.myelin.alldepths$depth <- substr(row.names(regional.myelin.alldepths), 1, 7) #assign depth variable
regional.myelin.alldepths$depth <- factor(regional.myelin.alldepths$depth, ordered = T, levels = ordered_depths)

regional.myelin.alldepths <-  regional.myelin.alldepths[!(regional.myelin.alldepths$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

regional.myelin.alldepths %>% 
  ggplot(aes(x = mean.R1, y = depth)) +
  geom_violin(aes(fill = depth), color = "white") +
  stat_summary(fun = median, geom = "point", shape = 23, size = 2.25, color = "white") +
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  scale_x_continuous(limits = c(0.495, 0.65), expand = c(0,0)) +
  scale_y_discrete(expand = c(0, 0)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depthplot_regionalR1.pdf", device = "pdf", dpi = 500, width = 2.6, height = 4.9)
for(this.depth in c(unique(regional.myelin.alldepths$depth))){
  depth.plot <- regional.myelin.alldepths %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp") %>%
  ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), colour=I("gray50"), size=I(.06)) + 
    paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.5, 0.62), oob = squish, na.value = "gray75") + theme(legend.position = "none")
  print(depth.plot)
  ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/%s_corticalmap.pdf", this.depth), device = "pdf", dpi = 500, width = 4.75, height = 2.25)
}

R1 Myelin Profile in Individual Regions

area4.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI")
area46.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI")
areaa24.myelin.alldepths <- regional.myelin.alldepths %>% filter(orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI")


ggplot() +
  geom_smooth(data = area4.myelin.alldepths, aes(x = mean.R1, y = depth, group = orig_parcelname), color = myelin.colorbar[4], se = F, size = .5) + 
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  scale_x_continuous(breaks = c(0.54, 0.57, 0.60, 0.63, 0.66), limits = c(0.54, 0.67)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 6, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 6, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area4_profile.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.25)


ggplot() +
  geom_smooth(data = area46.myelin.alldepths, aes(x = mean.R1, y = depth, group = orig_parcelname), color = myelin.colorbar[3], se = F, size = .5) + 
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  scale_x_continuous(breaks = c(0.51, 0.53, 0.55, 0.57, 0.59), limits = c(0.50, 0.605)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 6, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 6, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area46_profile.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.25)


ggplot() +
  geom_smooth(data = areaa24.myelin.alldepths, aes(x = mean.R1, y = depth, group = orig_parcelname), color = myelin.colorbar[2], se = F, size = .5) + 
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  scale_x_continuous(breaks = c(0.50, 0.52, 0.54, 0.56, 0.58), limits = c(0.495, 0.58)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Areaa24_profile.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.25)

Region-wise R1 Correlations Between Depths

regional.myelin.alldepths.wide <- regional.myelin.alldepths %>% pivot_wider(id_cols = orig_parcelname, names_from = depth, values_from = mean.R1) %>% select(-orig_parcelname)
regional.myelin.alldepths.wide <- regional.myelin.alldepths.wide %>% select(depth_7, depth_6, depth_5, depth_4, depth_3, depth_2, depth_1)
depth.R1.cors <- cor(regional.myelin.alldepths.wide) # compute correlation matrix

ggcorrplot(corr = depth.R1.cors, method = "square", type = "upper", show.diag = T, lab = F, lab_size = ) +
scale_fill_gradientn(colors = c("#edf6ff", "#c2d4ed", "#97adcc", "#345c91"))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depth_correlationmatrix.pdf", device = "pdf", dpi = 500, width = 4.4, height = 2.2)